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We study the problem of optimal control of dissipative quantum dynamics. Although under most 
circumstances dissipation leads to an increase in entropy (or a decrease in purity) of the system, 
there is an important class of problems for which dissipation with external control can decrease 
the entropy (or increase the purity) of the system. An important example is laser cooling. In 
such systems, there is an interplay of the Hamiltonian part of the dynamics, which is controllable 
and the dissipative part of the dynamics, which is uncontrollable. The strategy is to control the 
Hamiltonian portion of the evolution in such a way that the dissipation causes the purity of the 
system to increase rather than decrease. The goal of this paper is to find the strategy that leads to 
maximal purity at the final time. Under the assumption that Hamiltonian control is complete and 
arbitrarily fast, we provide a general framework by which to calculate optimal cooling strategies. 
These assumptions lead to a great simplification, in which the control problem can be reformulated 
in terms of the spectrum of eigenvalues of p, rather than p itself. By combining this formulation with 
the Hamilton- Jacobi-Bellman theorem we are able to obtain an equation for the globaly optimal 
cooling strategy in terms of the spectrum of the density matrix. For the three-level A system, we 
provide a complete analytic solution for the optimal cooling strategy. For this system it is found 
that the optimal strategy does not exploit system coherences and is a 'greedy' strategy, in which 
the purity is increased maximally at each instant. 
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I. INTRODUCTION 

In the last 15 years, optimal control theory (OCT) 
has been applied to an increasingly wide number of 
problems in physics and chemistry whose dynamics are 
governed by the time-dependent Schrodinger equation 
(TDSE). Theseproblems include control of chemical re- 
actions [1 HJl EL IIJEJI , state-to-state popula- 
tion transfer [aEfllllJll, 

shaped wavepackets [l3| . 
NMR s pin dynamics [lj], Bose-Einstein condensation 
[Hill 113, quantum computin gj lljl l |2C|. oriented 
rotational wavepackets [2l|. etc. |22l |23| . More recently, 
there has been vigorous effort in studying the control of 
systems governed by the Liouville-von Neumann (LVN) 
equation, where the central o bjec t is the den sity matrix, 
rather than the wavefunction pi El I2H I2I I2I l30| . 
The Liouville-von Neumann equation is an extension of 
the TDSE that allows for the inclusion of dissipative pro- 
cesses. Important examples of what may be thought of 
as quantum control processes that require the use of the 
LVN include laser control of chemical reactions in solu- 
tion, laser cooling, and coherence transfer in multi-spin 
systems. In all these cases, the external field (the laser or 
the RF field) is the coherent control, while the source of 
dissipation is contact with the environment. In the case 
of laser cooling, the environment is the vacuum modes of 
the electromagnetic field and the source of dissipation is 
spontaneous emission. 

In the majority of problems on control of quantum sys- 



tems dissipation is a nuisance; the purpose of the control 
is to either avoid, delay or cancel the dissipation pro- 
cess. Yet there is a remarkable exception to this pattern 
- laser cooling. The goal of laser cooling is expressed 
alternatively as increasing the phase space density, or de- 
creasing the entropy of the system. Purely Hamiltonian 
manipulations can in fact do neither, and therefore dis- 
sipation, rather than being a nuisance, is actually nec- 
essary to achieve true cooling. The optimal control of 
systems of this type is fascinating. The control itself, no 
matter what its time-dependence, leads only to Hamil- 
tonian evolution and hence no true progress toward the 
objective. On the other hand, the dissipation, while it 
is capable of producing progress toward the objective, is 
fundamentally not controllable and could in fact lead to 
a decrease in the objective. 

In ref. [2(|, we elucidated the interplay of the con- 
trolled, Hamiltonian evolution, and the uncontrolled, dis- 
sipative evolution in producing cooling. The "cooling 
laser" , while not directly cooling the system, in fact steers 
it to a region of parameter space where spontaneous emis- 
sion leads to cooling rather than heating. We define such 
a controlled manipulation as a "purity increasing trans- 
formation" . We believe that the study of such transfor- 
mations in their general mathematical context is of ex- 
treme interest, both in terms of discovering a wider class 
of physical processes where purity, and therefore coher- 
ence content can be increased, as well as because of the 
rich mathematical structure of the problems involving in- 
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terplay of Hamiltonian and dissipative dynamics. 

In [2|| , we solved the problem of optimal cooling for a 
2-level system completely, under the assumption of com- 
plete and rapid Hamiltonian control. We showed that the 
optimal cooling strategy in the 2-level system avoids pro- 
ducing coherences in the density matrix. Here we present 
a general framework for the analysis of optimal control 
in a system of N excited states coupled radiatively to M 
ground states, under the same assumptions. Using this 
framework we explicitly provide the optimal strategy for 
cooling of a three level A system. 

We first introduce the Lindblad dissipation model and 
a generalized concept of purity in section^ In Section 
III CI the problem of optimal cooling of a quantum me- 
chanical system is formulated. It is shown in Section IlIII 
that this problem can be reformulated solely in terms oi 
the eigenvalue distribution of the density operator. In 
doing this, we derive a reduced equation of motion for 
the spectral evolution under dissipation, parameterized 
by the unitary control (Sections IIII Bl and IIII CI) . Section 
II VI introduces the mathematical tools for finding optima] 
cooling strategies, namely the Hamilton- Jacobi-Bellman 
theorem. Section [V] provides an explicit description oi 
the optimal cooling strategy for the three level A sys- 
tem and proves its optimality. Finally we discuss future 
directions and conclude in Section IVll 



The second property follows from the fact that p = p* 
and therefore p{t) = p'(t). We will later derive an ex- 
plicit expression for the evolution of the spectrum of the 
density operator under dissipation. The third property 
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FIG. 1: A general N + M level quantum system with spon- 
taneous emission rates between various energy levels. 



II. SETTING UP THE CONTROL PROBLEM 

A. The system equations of motion and the 
Lindblad formula for dissipation 

Let p denote the density matrix of an N level quantum 
system (see figure The density matrix evolves under 
the Liouville von Neumann (LVN) equation which takes 
the form 

p = -i[H{t),p]+L{p) (1) 

where — i[H, p] is the unitary evolution of the quantum 
system and L(p) is the dissipative part of the evolution. 
The term L(p) is linear in p and is given by the Lindblad 
form 0,113, i.e. 

ij 

where Fij are the Lindblad operators. In this manuscript, 
we assume the only relaxation mechanism is spontaneous 
emission and therefore we take Fij = ^jij Eij where the 
operator Eij = and 7y represents the rate of spon- 
taneous emission from level j to level i. Eq. Q has the 
following three well known properties: 1) Tr(p) remains 
unity for all time, 2) p remains a Hermitian matrix, and 
3) p stays positive semi-definite, i.e. that p never devel- 
ops non-negative eigenvalues. 
The first property follows from 



B. Definitions of Purity 

The density matrix is capable of describing any mixed 
state in quantum mechanics, ranging from pure states 
that are solutions of the TDSE, to completely incoherent 
states. There are several common ways of characterizing 
how close an arbitrary mixed state p is to a pure state. 
These measures can be generally termed purity measures 
or purities. We use P{p) to denote the purity of the 
density operator p. 

The most common, and perhaps the simplest measure 
is Tr(p 2 ) [HIH HIIH. For any density matrix, < 
Tr(p 2 ) < 1, with equality only for a pure state. Thus, 
the larger the value of Tr(p 2 ), the closer a state is to 
being pure. Another useful measure is the von Neumann 
entropy, Svn — ~k Tr(plnp) j^. The von Neumann 
entropy goes to zero for a pure state and is greater than 
zero for any mixed state, and thus the size of the von 
Neumann entropy is a measure of the degree of impurity 
of a state. Two other measures are the largest eigenvalue 
of p, \p\oo, which goes to 1 for a pure state and is less 
than one for a mixed state; and a measure based on the 
expansion of the characteristic equation for p, which has 
Tr(p 2 ) as its leading term, but also takes into account 
higher order terms, e.g. Tr(/3 3 ) l^. 1 



Tr(p) = Tr(-i[H,p]) + Tr(Lp) = 0. (2) 



The purity function can be thought of mathematically as a 
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In general, as is apparent from the above discussion, 
the entire density matrix p is not needed in order to char- 
acterize the purity of the system; rather, all that is nec- 
essary is the set of eigenvalues A of p. All purities can 
therefore be defined as functions solely of the eigenvalues, 
i.e 

P(p) = P(X(p)). (4) 

We will use the following definition of purity for the 
remaining part of the paper. 

Definition 1 Given the density operator p, with spec- 
trum X, define its purity P(p) as the largest eigenvalue 
of p, i.e. 

P{ P ) = \p\ 00 = lim Tr{p n ) * = X[. (5) 

n — >oo 

Here A^ is the vector of eigenvalues of p arranged in a 
decreasing order; for the remainder of this paper the su- 
perscript I will be assumed every time A is written. Al- 
though many of the results in this paper are very general, 
we choose this measure as it gives simple answers for the 
cooling strategies. We will often use P(p) or P(X) to 
mean the same thing, where it is understood that A is 
the spectrum corresponding to p. 

C. Formulation of the Control Problem 

The problem we address in this paper is the con- 
trol of purity content of a quantum dissipative system 
which evolves under the LVN equation of motion given by 
eq. The Hamiltonian H[E] depends on an externaly 
controlled laser field E(t) through the dipole coupling 



(partial) ordering over the set of allowed eigenvalues such that 
the totally pure state having the spectrum [1,0, ...,0] yields the 
greatest value of purity and the totally mixed state with spec- 
trum [i i, i] yields the lowest. A necessary minimum 
of structure on t he puri ty ordering is provided by the concept 
of majorization |37t I'tti . Let x and y be two d-dimensional 
real vectors. We use the notation to indicate the vector 
whose entries are the entries of x, arranged into decreasing order, 
Xj > > ■ ■ • > xjj. We say x is majorized by y, written x -i y, 
if 

V-'V N>;. ( 3 ) 

3=1 3=1 

for k = 1, • • ■ , d, with equality when k = d. Loosely speaking, 
this definition gives quantitative meaning to the amount of dis- 
order or mixing in a collection of real numbers. For example, 
for any < t < 1, [^, ^] -< [t, 1 — t]. For any d dimensional 
probability distribution p, [i, ■ ■ • , ^] -i [pi, • • ■ , p<j]. Note that 
there are vectors x and y which are incomparable in the sense 
that neither x -< y nor y -< x (for example x = [0.5,0.25,0.25] 
and y = [0.4,0.4,0.2]); majorization therefore gives only partial 
ordering. Any reasonable measure of purity should respect the 
majorization relation, namely for two eigenvalue distributions 
we should have P(A') < P(A") if A' -< A". Such functions are 
termed Schur-convex. 



term. Beginning with the system in an initial mixed 
state it is required to find a control field functionality 
E{t) that will drive the system through its equations of 
motion to maximal purity, as defined by eq. JSJ), at 
some final time T . 

The system evolution equation contains both a Hamil- 
tonian part, 

p H = -i[H[E],p], 
and a dissipative part, given by 
Pd = L(p). 

The Hamiltonian term leads to unitary evolution which 
does not change the spectrum, and the purity depends 
only on the spectrum. Thus, the dissipative term is re- 
quired to obtain a purity increase. In pfij . the control 
problem was solved completely for the two-level system. 
In this paper we develop a formalism applicable to gen- 
eral A^-level systems. 

III. REFORMULATION OF THE CONTROL 
PROBLEM IN TERMS OF THE SPECTRUM OF p 

A. Simplifying assumptions: complete and 
instantaneous unitary control 

In this section we develop a general formalism that 
highlights the cooperative interplay between Hamiltonian 
and dissipative dynamics. Following |26| . we assume that 
the action of the control Hamiltonian can be produced on 
a time scale fast compared with spontaneous emission. 
This assumption is well established on physical grounds, 
since femtosecond laser control is now widely available 
and typical spontaneous emission times are nanosecond. 
In this paper we make an additional useful simplifying 
assumption about the dynamics, namely that the control 
Hamiltonian H{t) can produce any unitary transforma- 
tion U S SU (N) in the N level system, i.e. the system of 
interest is unitarily controllable. Combining these two as- 
sumptions we have that any unitary transformation can 
be produced on the system in negligible time compared 
to the dissipation. 

We use the notation 

0(p) = {UpU*\u eSU(N)}, 

to denote the orbit of p under unitary transformations. 
Since X(UpW) = X(p), it is obvious that P(p) = P(X) 
is constant along the orbit 0(p); however P is not: the 
rate of change of the purity due to dissipation is affected 
by where in 0(p) the density matrix resides. In other 
words, due to the 'instantaneous controllability' assump- 
tion, unitary controls can instantaneously direct p along 
the orbit in order to change P in a controlled manner. 

The above dynamical assumptions lead to another very 
important simplification. Since we have assumed that 
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all unitary transformations in SU(N) can be produced 
instantaneously, this includes bringing the density matrix 
into diagonal form. As a result, the different elements of 
each orbit can be considered redundant, and the orbit of 
p can be completely represented by its diagonal form, or 
'spectrum', A(p). This suggests reformulating the control 
problem entirely in terms of the spectrum, rather than 
in terms of p itself. The key step in this reformulation 
is to replace the equation of motion for p, eq. QJ, with 
an equation of motion for the spectrum. We do this in 
the next section. As the purity is a function solely of the 
spectrum, this equation will allow the optimization to be 
performed just on the set of allowed spectra, significantly 
reducing the complexity of the problem. The controls 
will enter into this equation in a modified way that gives 
additional insight into the interplay of Hamiltonian and 
dissipative dynamics. 

B. Equations of Motion for the Eigenvalues 
Assuming Fast Unitary Evolution 

Suppose that p has a nondegenerate spectrum, and let 
A be its associated diagonal form. Consider two unitary 
transformations, U\ and U^- Then both p\ — V\AJJ\ 
and P2 — U2AU2 belong to 0(p). However, they do not 
have the same spectrum after evolution under the dis- 
sipative dynamics. To understand how the spectrum of 
the density operator evolves, note that Hamiltonian dy- 
namics produces no change in the spectrum. Therefore, 
the change in the spectrum is solely due to dissipation. 
After small time St the initial density operator p evolves 
to 

p^p + L(p)St, (6) 

If A represents the diagonalization of the original den- 
sity operator p (p = UAW , where U is unitary) then the 
new density operator can be written as 

p -> p + L(p)5t = U(A + tfL(UMJ*)U8t)U\ (7) 

Consider now the change in spectrum under the evolu- 
tion of eq. iJJJl ■ Since A is diagonal, the spectrum on the 
right hand side is, to first order in St, just the diagonal 2 
i.e. 

X(t + St) = diag(A + [7 f L(UAU r )USt). 

Given the matrix A, the notation diag(A) represents a 
vector whose entries are the diagonal entries of A. The 



2 This is simply the well known result of first order perturbation 
theory which, when applied to a perturbed Hamiltonian, states 
that the first order corrections to the energies are the diagonal 
elements of the perturbing Hamiltonian V. 



rate of change of eigenvalues is then 

A = diag{U ] ' L{U AU f )U) (8) 

which is in general different for different choices of U. 
Thus by applying varying unitary transformations U and 
letting the dissipative dynamics evolve for some small 
time St we get different evolution of the spectrum. The 
unitary transformation should therefore be thought of as 
a control by which the spectrum of the density matrix 
can be affected. 



C. Canonical decomposition 

To proceed further, observe that the right hand side 
of eq. (JHJ) describing the change in the spectrum under 
operation of the Lindbladian is a linear transformation 
on the vector of eigenvalues (see appendix |A"|) 

A = MA. 

To obtain an explicit expression for M first note that for 
U = I in eq. (JHJ we have A = AX with A a Q-matrix 
(columns sum to zero) defined by Aij — 7^ for i ^ j and 
Mi — - J2k 7k» otherwise. We split 

A = B + D 

where D is the diagonal part of A and is all non-positive 
whereas B contains all off-diagonal entries and is all non- 
negative. Using these definitions we get for general U in 
eq. (JHJ (for details see appendix |A"|) : 

A = (e T B9 + e T oD)A, (9) 

where @ij = \Uij\ 2 , is the Schur product of U with 
its complex conjugate. Note that O has the important 
property of being a doubly-stochastic matrix (rows and 
columns all sum to unity). The notation T o D de- 
notes the linear transformation of the diagonal of D (as 
a vector) under the action of T . In other words, if d = 
diag(D), then Q T o D is a diagonal matrix whose diago- 
nal is Q T d. Note that in the special case where U £ {Pi} 
— the set of permutations — = Pi, Pf o D = Pf DP L 
and hence eq.© simplifies to A = T /10A. 

Eq. © is one of the central results of this paper; it 
provides a reduced equation of motion for the spectral 
evolution under Lindblad dissipation and parametrized 
by the unitary control. From eq. ©, it is straightforward 
to infer, for example, that the eigenvalues of the density 
operator always remain nonncgative. In order to become 
negative an eigenvalue must pass through zero. If any of 
the eigenvalues \j — 0, however, the only contributions 
to Xj will be nonnegative since the only nonpositive ele- 
ments in M — Q T BQ + Q T o D reside on the diagonal. 
Hence none of the eigenvalues can turn negative. 
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D. Revised definition of the Control problem 

Having formulated an equation of motion for the spec- 
trum, eq. © , we can now redefine the control problem in 
terms of the spectrum alone. We seek a control strategy 
in the form of a time varying unitary-stochastic matrix 
O(i) which when applied to the spectral equation of mo- 
tion 0, will produce maximal purity P(X) at the final 
time T. 

One strategy for choosing Q(t) is to instantaneously 
maximize the purity P(X) at each point in time. Maxi- 
mization algorithms that utilize this strategy are termed 
'greedy' algorithms and do not in general guarantee ob- 
taining maximum possible purity at the final time T. To 
calculate the globally optimal coolin g st rategy we use the 
principle of dynamic programming |39| . as described in 
the next section. 



IV. DYNAMIC PROGRAMMING AND THE 
HAMILTON- J ACOBI-BELLMAN PDE 

We now use the principle of dynamic programming for 
finding the optimal Q(t) in Eq. 0. We will develop the 
basic ideas through the problem under consideration. Let 
V(A, t) denote the maximum achievable purity starting 
from initial eigenvalue spectrum A at time t (T — t units 
of time remaining). By definition of V(X,t), it is the 
maximum achievable purity if O is chosen optimally over 
the interval [t, T]. Suppose that at time t, the spectrum 
of p(t) is X(t) and we make a choice of 6(t) . The resulting 
density operator after time St depends on the choice of 
0(f). The choice of 0(f) should be such that for the 
resulting new spectrum A(t + St), the return function, 
V(X(t + St),t + St) is maximized and by definition of the 
optimal return function should be same as V(X(t), t). By 
a Taylor series expansion we obtain 

V(X(t + St),t + St) = v(x(t),t) + dv { x,t) st 

ot 

+ max(^M £A(0)).(1O) 

This then gives the well known Hamilton- Jacobi-Bellman 
PDE 

dV(X,t)_8V(X,t) , w ^ )A(e))=( , (11) 



dt 



Ot 



e dX 



Observe that at the final time T, the value of the return 
function is just the purity of the density operator, i.e. 

V(X,T)=P(X). 

If we solve this PDE, together with its final condition, 
we will get the optimal control as a function of the 
spectrum A and the time t, denoted as = 0*(A, t). 
In other words, given the spectrum A of the density op- 
erator at time t, the best cooling strategy is to choose 
<d*(X,t). This implies that the control problem is solved 



not just for a particular set of initial conditions; rather, it 
is embedded in a wider problem and a solution is sought 
simultaneously for all possible initial conditions. 
In equation (|llf> . the term SV gt'^ has no dependence on 
0, therefore 

0*(A,i) =argmax(^M A(O)). 
Substituting for A(0) from Eq. ©, yields 



0*(A, t) = argmax 



dV(X ' t \(Q T BQ + QoD)X). (12) 



e v dX 

Thus the problem reduces to finding the optimal control 
0*(A,i) that maximizes the expression 



F(0) = p T (Q T B& + & 1 o D) X. 



(13) 



where the vector p, is defined as pj — ^j- (Although p is a 
function of A and t, we just use p and keep in mind that 
the dependence is implied). Note that a priori V(X,t) 
and hence p are not known. However if we can make 
a guess at the optimal control strategy (which depends 
on A and t) and use this optimal strategy to integrate 
the equation of motion of the system evolution to obtain 
V (A, t) and hence p, then we can verify if the optimal 
control and the corresponding p satisfy equation i|12fl . 
We illustrate this by finding optimal cooling strategies 
for a 3-level Lambda system. 

The following properties of equation ljl3|l will be used 
subsequently. being a double stochastic matrix implies 
that [1, 1, ...1]0 = [1, 1, ...1]. Furthermore [1, 1, ...1](B + 
D) = and therefore F(&) vanishes for p T = [1, 1, 1]. 
The elements of p can therefore be shifted by a constant 
amount to make a specific component of p vanish without 
influencing the value of F(0). 



V. SOLUTION OF THE OPTIMAL CONTROL 
PROBLEM FOR THE 3 LEVEL SYSTEM 

A. Preliminaries 

Consider a three level Lambda system depicted in fig- 
ure |3 The excited state spontaneously decays into the 
stable ground states |1) and |3) at rate 71 and 72 respec- 
tively. We will assume without loss of generality that 
7i > 72- 

The evolution of the density matrix of the three level 
Lambda system is given by 

p = -i[H(t),p] + <n(E lP El - ±{ElE 1)P }) 

+ l2 {E 2 pEl- l -{E\E 2 ,p}), (14) 

where E x = |1)(2| and E 2 = |3)(2|. The equation of 
motion for the spectrum of the density matrix is then @ 
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FIG. 2: A 3 level Lambda system with spontaneous emission 
rate from level 2 to 1 given by 71 and spontaneous emission 
rate from 2 to 3 given by 72. 
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FIG. 3: The eigenvalue evolution under the optimal cooling 
strategy for the three level A system. 
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The objective is to maximize the purity at time T, P(T), 
as measured by the largest eigenvalue of p (Definition 1). 



B. The optimal strategy: Keep p diagonal and 
ordered 

Given the equation of motion defined by eq. ^]and the 
objective defined by Definition 1, we have the following 
theorem: 

Theorem 1 The optimal cooling strategy For the 

3-level system described above (labelled as in fig. EJ) 
if any unitary transformation U G SU (3) can be pro- 
duced in arbitrarily small time, then the optimal cool- 
ing strategy is to keep the density operator p(t) diagonal 
for all times (produce no coherences) and ordered i.e. 
Pn(t) > p 22 (t) > p 33 (t). 

The optimal control strategy has the following alternate 
description. Throughout the cooling process, we keep the 
largest eigenvalue in the eigenstate |1), the next largest 
in state |2) and finally the smallest in state |3). As the 
population in state \2) decays spontaneously to state |1) 
and 1 3), after some time t* , the population of states |2) 
and 1 3) will become equal. From that point onwards, 
we always maintain the population of states |2) and |3), 
equal (see figure |3J). We will refer to this strategy as 



"greedy" since it maximizes the rate of increase of the 
objective at each point in time. 

To prove optimality of the above strategy we proceed 
as follows. We first compute V(X,t) for the proposed 
strategy and then show that it satisfies the HJB equa- 
tion maximized over all unitary transformations. Follow- 
ing the convention that the elements of the vector A are 
arranged in decreasing order, this amounts to showing 
that 

/ = argmax F(Q) 

ee{\u\ 2 l uesu(N)} 

■ T (e T Be + e T od) A{i5) 



argmax p 
ee{\u\ 2 l u<£SU{N)} 



where / is the identity operator. This implies that the 
eigenvalues should be continuously maintained in their 
ordered arrangement. Note that despite the simplicity 
of this result, in general the continuous intervention of 
a control field is required in order that this condition be 
fulfilled. 



C. The Return Function for the Ordered Diagonal 
Strategy 

We now evaluate the return function for the putative 
optimal strategy. Let r = T—t denote the remaining time 
for cooling. According to the strategy proposed above, 
two evolution regimes exist depending on whether t < t* 
or t > t* , where r* is the critical time required for A2 
and A3 to come to equilibrium. 

In the case where r < t* , under the proposed strategy 
the evolution equations of the system take the form 



Ai = 71A2 

A 2 = -(71 + 72) A 2 

A3 = 72 A 2 



(16) 
(17) 
(18) 
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and therefore 
V(X,T-t) = \ 1 



] A 2 {l- e - { ^ 2)T \ : T ■:. T 



7l + 72 



(19) 

By definition A 2 (r*) = A 3 (r*) and A 2 (r*) = 
/ \ 2 g-(7i+72)r _ Using these equalities, the following ex- 
plicit forms for A 2 (r*) and r* can be computed 



A 2 (r*) 



72A2 + (71 + 72)A 3 
7i + 272 



7i+72 V A 2 (7i + 27 2 ) J 

After this point in time, under the ordered diagonal pol- 
icy, the populations of states |2) and |3) are maintained 
at equilibrium such that A2(r) = A 3 (r) = |(1 — Ai(r)). 
The system dynamics therefore takes the form 

from which the return function for the regime r > r*, 
can be explicitly computed. 



V(X, T —t) = 1 - 2A 2 (r*)e~- 



T > T 



(21) 



As the return function enters the HJB equations only 
through its derivatives fi — 4|j, we proceed to compute 
these derivatives explicitly for use in the next section. 
For t < r*. we have 



Mi = 1 



M2 = 



7i 



_ g-(7l+72)r^ 



(22) 



7i + 72 
M3 - 
and for r > t* , we have 
Mi = 

272A2 + 71 A 3 Q( T _ T ») 
A 2 (7i + 272) 



M2 



M3 



Note that in both regimes fii > fi 2 > M3 1 a property that 
will be used below. 3 Also note that the /x's are continuous 
at t — t* up to a constant shift (see remark at the end 
of section HV|) . 



D. Proof that the Return Function for the Ordered 
Diagonal Strategy Satisfies HJB 

We proceed to calculate argmaxF(O) for F(&) given 
by eq. (JT3J and // given by eq. and We show 

that 9* = 7 and hence the ordered diagonal strategy 
satisfies the HJB equation, proving that this strategy is 
globally optimal. 



3 In order to prove this statement for r < t* note 
that in this regime A3(r) < \2(t), which implies 

A 2 



a. Absence of ground state coherences in the ordered 
diagonal solution We first prove that the optimal trans- 
formation in equation (jlHfl has the property that 
©13 = ©3i = 0, namely that the ground state coherences 
vanish throughout the evolution of the optimal trajec- 
tory. Suppose 0i3 7^ and ©31 ^ and say ©31 > ©13. 
From equation l|13|l we have 

F(Q) = [7i(mi©11 + M3©13) + 72(Ml©31 + M3©33)] 
X [Ai©21 + A 2 ©22 + A 3 ©23] 
- (71 + 72)[MlAl@21 + M3A 3 ©23], 

(24) 

where we have chosen \i 2 = and hence [i\ > > /.*3. 
Let A = 13 . Observe that in the above equation we 
can increase ©n and O33 by an amount A and decrease 
13 and ©3i by A, to generate a new doubly stochastic 
matrix which gives a larger value of F(Q) (this follows 
from the relations 71 > 72 and \i\ > > fj,a). Hence we 
assume ©i 3 = 0. Let Ai be the new value of ©31. Now 
if we increase ©n and ©32 by Ai and decrease 631 and 
©12 by the same amount we get a new doubly stochastic 
matrix which gives a larger value of F(Q). Hence we 
need to maximize ^(0) only over those doubly stochastic 
matrices for which ©13 = 3 i = 0. 

b. Dependence of F(Q) on the remaining parameters 
in O As the rows and columns of must sum to unity 
there remain only two degrees of freedom in the compo- 
nents of ©. Therefore, we can write F(Q) as a function 
of only two of its components: 

F(0) = F(0 21 ,0 23 ) 

= [71/11(1 - ©21) +72M3(1 - ©23)] 
x [A 2 + 2 i(Ai - A 2 ) + © 23 (A 3 - A 2 )] 

-(7l +72) [Ml@2lAl +M3@23A 3 ] ■ (25) 

It is now required to find the maximum of ^(©21, ©23) 
on the triangular domain < ©21 < 1, < ©23 < 1, 

©21 +©23 < 1. 

c. The maximum cannot lie at an interior point. 
Suppose F has a maximum in the interior, then the Hes- 
sian of F at that point must be negative definite. We 



proceed to show that the Hessian Gi 



with 



71+72 



((71 +2 72 )e-(71+72)r _ 72 ) > X 3 . 



i, j = {1,3}, is not negative definite anywhere and there- 
fore the maximum must reside on the boundary. Com- 
puting the components of G we find 

Gn = -2(Ai - A 2 )mi7i 
G33 = — 2(A 3 — A 2 )^372 

G31 = G13 = -Mi(A3 - A 2 )7i - M3(Ai - A 2 )7 2 -(26) 

Denoting a = 71^1 (A3 — A2) and b = 72M3(Ai — A2), the 
determinant of G is 4ab — (a + b) 2 = — (a — b) 2 < 
such that one of the eigenvalues of G is non negative and 
therefore G is not negative definite. 

d. The maximum point is (©21,623) = (0,0). As 
the maximum does not reside in the interior of the trian- 
gular domain it must lie on one of the edges [(0,0),(0, 1)], 
[(0,0), (1,0)] or [(0,1), (1,0)]. 
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It can be shown by checking the first and second deriva- 
tives along the edge [(0, 1), (1, 0)] that the maximum 
along that interval lies at the end point (021,823) = 
(0, 1). We now check the remaining two edges. As Gn 
and G33 are both non positive it follows that F(Q) is 
concave in both the O21 and 023 directions. Therefore, 
if in addition the slope at the point (0, 0) is negative in 
both directions, this establishes the existence of a maxi- 
mum at that point. We proceed to show that indeed the 
slopes are non negative: 



dF 



8Q 21 
dF 



23 



(0:0) 



(0.0) 



(Ai - A 2 )[^372 + Mi7i] + W7i] 
-(71 +72V1A1 < (27) 
(A 3 - A 2 )[^ 3 72 + Mi7i] + M-M373] 
-(71 +72V3A3 < 0. (28) 



The first expression follows from the fact that 71 > 72, 
/ii > > /U3 and Ai > A2. The second expression can be 
proved by inserting the explicit forms for /z, eq. I|22|) and 
(|23[1 . for the two regimes of T — t. 



VI. DISCUSSION AND CONCLUSIONS 

We have presented a general framework for calculating 
optimal purity increasing strategies in N level dissipative 
systems under the assumption of complete and instanta- 
neous unitary control. In so doing, we derived a reduced 
equation of motion for the spectral evolution under dis- 
sipation and parametrized by the unitary control. The 
Hamilton- Jacobi-Bellman Theorem was invoked to pro- 
vide sufficient criteria for global optimality. This general 
framework was then explicitly applied to derive and prove 
optimality of the greedy cooling strategy for a three level 
A system. 

In future work we intend to apply this methodology to 
obtain explicit optimal cooling strategies for general N + 
M level systems comprised of M excited states coupled to 
iV ground states. One is tempted, by extrapolation from 
the present results, to assume that the greedy algorithm 
should be optimal in general and hence that coherences 
do not play a role in the optimal cooling strategy for 
N > 3. However, preliminary numerical results based on 
dynamical programming show that the greedy algorithm 



is in general not optimal in these systems. Rather, a 
strategy based on "delayed gratification" is superior to 
the greedy strategy, and coherences play a small but finite 
role in these larger systems. This will be the subject of 
a future paper. 

APPENDIX A: DERIVATION OF THE 
'CANONICAL FORM' 

We wish to show that eq. (jHJ is a linear transformation 
of the form 

A = MA. 

Recall first that 
U^L(UAU r )U 

= E^ C/t \\W\ UA = tfm\ - \ {\j)(j\,UAU*} I 



tf\i){j\UAtf\j){i\U - -{tf\j){j\U,A} 



Expanding eq. JSJ and rewriting it in component form 
we have 



s 



E 



^M ks \ s , 



(AI) 



with 



m = e T Be + e 1 o d, 



and with the definitions of O, B, D and the operation 
O o D as provided in the main text. Rewriting the above 
in vector format we have precisely eq. @ 

A = {Q T BQ + e T oD)X. 
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